Investigation of the apical progenitors diversity

Load libraries and QCFiltered dataset

# Load library
library(Seurat)
library(monocle)
library(princurve)
library(cluster)
library(parallel)
library(scrattch.hicat)
library(scrattch.vis)
library(ggplot2)
library(ggExtra)
library(dplyr)
library(reshape)
library(gridExtra)
library(RColorBrewer)
library(wesanderson)
library(viridis)

#Set ggplot theme as classic
theme_set(theme_classic())

# Load the full annotated dataset
Allcells.data <- readRDS("./QC.filtered.cells.RDS")

Extract Apical progenitors and compute Pseudo-DorsoVentral score

Extract apical progenitors and filter the gene expression matrix

We perform Kmeans clustering on the 2 cell state scores :

• Apical progenitors AP
• Basal progenitors BP

set.seed(100)
#K-means clustering based on AP, BP scores across cells
cl <- kmeans(cbind(Allcells.data@meta.data$AP_signature1, Allcells.data@meta.data$BP_signature1), 3)
Allcells.data@meta.data$kmeanClust <- paste0("Clust.",cl$cluster)

col.pal <- wes_palette("GrandBudapest1", 3, type = "discrete")

p1 <- ggplot(Allcells.data@meta.data, aes(x=AP_signature1, y=BP_signature1, colour = kmeanClust)) +
  scale_color_manual(values=col.pal) +
  geom_point() + 
  theme(legend.position="none")
ggMarginal(p1, type = "histogram", fill="lightgrey") ; rm(p1)

DimPlot(Allcells.data,
        group.by = "kmeanClust",
        reduction.use = "spring",
        cols.use = col.pal,
        dim.1 = 1,
        dim.2 = 2,
        do.label=T,
        label.size = 4,
        no.legend = F)

We then extract the glutamatergic neuron branche as beeing the Kmeans cluster with the highest mean Apical progenitor signature

#Fin cluster wiht the highest mean APscore
MeanKclust.APscore <- aggregate(AP_signature1 ~ kmeanClust, Allcells.data@meta.data, mean)
APclust <- MeanKclust.APscore %>% filter(AP_signature1 == max(AP_signature1)) %>% pull(kmeanClust)

#Extract apical progenitors 
barcodes <- Allcells.data@meta.data %>% filter(kmeanClust == APclust) %>% pull(Barcodes)
AP.data <-  SubsetData(Allcells.data, cells.use = barcodes , subset.raw = T,  do.clean = F)

#Further filter the 3 outlier cells based on spring coordinates
cells <- rownames(AP.data@dr$spring@cell.embeddings[AP.data@dr$spring@cell.embeddings[,2] > 250,])
AP.data <-  SubsetData(AP.data, cells.use = cells , subset.raw = T,  do.clean = F)

DimPlot(AP.data,
        group.by = "kmeanClust",
        reduction.use = "spring",
        cols.use = col.pal,
        dim.1 = 1,
        dim.2 = 2,
        do.label=T,
        label.size = 4,
        no.legend = F)

Filter gene counts matrix

#Remove non epressed genes
num.cells <- Matrix::rowSums(AP.data@data > 0)
genes.use <- names(x = num.cells[which(x = num.cells >= 20)])
AP.data@raw.data <- AP.data@raw.data[genes.use, ]
AP.data@data <- AP.data@data[genes.use, ]

#Normalize and Scale the data
AP.data <- NormalizeData(object = AP.data,
                         normalization.method = "LogNormalize", 
                         scale.factor = round(median(AP.data@meta.data$nUMI)),
                         display.progress = F)

AP.data <- FindVariableGenes(object = AP.data,
                                  mean.function = ExpMean,
                                  dispersion.function = LogVMR,
                                  x.low.cutoff = 0.0125,
                                  x.high.cutoff = 3,
                                  y.cutoff = 1, do.plot = F,
                                  display.progress = F)

AP.data <- ScaleData(object = AP.data, vars.to.regress = c("CC.Difference","percent.mito", "nUMI"), display.progress = F)

Import Spring dimensionality reduction

Cell cycle associated genes were excluded for PCA dimensionality reduction and Spring plot was generating with these parameter :

Number of cells: 1648
Number of genes that passed filter: 857
Min expressing cells (gene filtering): 3
Min number of UMIs (gene filtering): 3
Gene variability %ile (gene filtering): 90
Number of principal components: 7
Number of nearest neighbors: 20
Number of force layout iterations: 500

#Import Spring coordinates calculated without cell cycle genes
Coordinates <- read.table("./Progenitors/E12.AP.Coordinates.txt", sep=",", header = F)[,c(2,3)]
rownames(Coordinates) <- rownames(AP.data@meta.data)
AP.data <- SetDimReduction(AP.data,
                           reduction.type = "spring.AP",
                           slot = "cell.embeddings",
                           new.data = as.matrix(Coordinates))
AP.data@dr$spring.AP@key <- "spring.AP"
colnames(AP.data@dr$spring.AP@cell.embeddings) <- paste0(GetDimReduction(object = AP.data, reduction.type = "spring.AP",slot = "key"), c(1,2))

Fit a principal curve over the AP in the Spring space

data <- data.frame(springAP.1 = AP.data@dr$spring.AP@cell.embeddings[,1],
                   springAP.2 = AP.data@dr$spring.AP@cell.embeddings[,2],
                   spring1 = AP.data@dr$spring@cell.embeddings[,1],
                   spring2 = AP.data@dr$spring@cell.embeddings[,2])

# Fit the principal curve
fit <- principal_curve(as.matrix(data[,1:2]),
                       smoother='lowess',
                       trace=TRUE,
                       f = 0.7,
                       stretch=0,
                       plot_iterations = F)

## Starting curve---distance^2: 2022449778
## Iteration 1---distance^2: 745172
## Iteration 2---distance^2: 696270.5
## Iteration 3---distance^2: 677509.1
## Iteration 4---distance^2: 668251.7
## Iteration 5---distance^2: 663393.6
## Iteration 6---distance^2: 660928.6
## Iteration 7---distance^2: 659586.4
## Iteration 8---distance^2: 658917.8
## Iteration 9---distance^2: 658638.8

DorsoVentral.Score <- fit$lambda/max(fit$lambda) #The actual speudotime
pc.line <- as.data.frame(fit$s[order(fit$lambda),]) #The principal curve smoothed

data$Phase <- as.character(AP.data@meta.data$Phase)
data$DorsoVentral.Score <- DorsoVentral.Score

# Direction of the maturation score using Zbtb20 expression (reverte if positive correlation)
if (cor(data$DorsoVentral.Score, AP.data@data['Zbtb20', ]) > 0) { data$DorsoVentral.Score <- -(data$DorsoVentral.Score - max(data$DorsoVentral.Score))}

AP.data@meta.data$DorsoVentral.Score <- data$DorsoVentral.Score

#Plot Cell onto PC1 and PC2 with principal curve
ggplot(data, aes(springAP.1, springAP.2)) + 
  geom_point(aes(color=Phase), size=2, shape=16, values=col.pal) +
  geom_line(data=pc.line, color='red', size=0.77)

#Plot Speudotime color gradient on the cell cycle filtered Spring embbeding
ggplot(data, aes(springAP.1, springAP.2)) +
  geom_point(aes(color=DorsoVentral.Score), size=2, shape=16) +
  scale_color_viridis(direction = -1) +
  geom_line(data=pc.line, color='red', size=0.77)

Manuscript Fig. 5A

#Plot Speudotime color gradient on the Spring embbeding calculated from full dataset
ggplot(data, aes(spring1, spring2)) + 
  geom_point(aes(color=DorsoVentral.Score), size=2, shape=16) +
  scale_color_viridis(direction = -1) 

Find differentially expressed genes along the pseudo DV axis

Initialize a monocle object

# Transfert metadata 
meta.data <- data.frame(barcode = rownames(AP.data@meta.data),
                        Cluster = AP.data@meta.data$old.ident,
                        DorsoVentral.Score =  AP.data@meta.data$DorsoVentral.Score,
                        CellcyclePhase = AP.data@meta.data$Phase,
                        row.names = rownames(AP.data@meta.data))
                   
Annot.data  <- new('AnnotatedDataFrame', data = meta.data)

# Transfert count data
count.data = data.frame(gene_short_name = rownames(AP.data@raw.data),
                  row.names = rownames(AP.data@raw.data))

feature.data <- new('AnnotatedDataFrame', data = count.data)

# Create the CellDataSet object
gbm_cds <- newCellDataSet(as.matrix(AP.data@raw.data),
                          phenoData = Annot.data,
                          featureData = feature.data,
                          lowerDetectionLimit = 1,
                          expressionFamily = negbinomial())

gbm_cds <- estimateSizeFactors(gbm_cds)
gbm_cds <- estimateDispersions(gbm_cds)
gbm_cds <- detectGenes(gbm_cds, min_expr = 0.1)

rm(list = ls()[!ls() %in% c("AP.data", "gbm_cds")])

Test each gene trend over pseudo-DV score

# Exclude cell cycle associated genes
CCgenes <- as.character(read.table("./Progenitors/CellCycleGenes.csv", sep = "\t", header = F)[,1])
Input.genes <- AP.data@var.genes[!AP.data@var.genes %in% CCgenes]

# Perform the test for differential expression as a function of pseudo-DV score while controling for cell cycle phase
DV.Axis.genes <- differentialGeneTest(gbm_cds[Input.genes,], 
                                      fullModelFormulaStr = "~sm.ns(DorsoVentral.Score, df = 3)*CellcyclePhase", 
                                      reducedModelFormulaStr = "~CellcyclePhase", 
                                      cores = detectCores() -2)

# Filter genes with a FDR < 0.001
DV.Axis.genes.FDR.filtered <- DV.Axis.genes %>% filter(qval < 1e-3)

Smooth significative gene expression pseudo-DV axis

# Create a new pseudo-DV vector of 500 points
nPoints <- 500
new_data <- data.frame(DorsoVentral.Score = seq(min(pData(gbm_cds)$DorsoVentral.Score), max(pData(gbm_cds)$DorsoVentral.Score), length.out = nPoints))

# Smooth gene expression
Smooth.curve.matrix <- genSmoothCurves(gbm_cds[as.character(DV.Axis.genes.FDR.filtered$gene_short_name),],
                                       trend_formula = "~sm.ns(DorsoVentral.Score, df = 3)",
                                       relative_expr = TRUE,
                                       new_data = new_data,
                                       cores= detectCores() - 2)

Cluster cell by similar expression profiles on the pseudo-DV axis

set.seed(100)
# Cluster cells using the Partitioning Around Medoids algorithm
Cells.Clust <- pam(as.dist((1 - cor(Smooth.curve.matrix ,method = "spearman"))), k=7)
Domaines.Clust <- data.frame(Domaines = paste0("Clust.",Cells.Clust$clustering))

Cluster gene by similar profiles over pseudo-DV axis

Perform hierarchical clustering

# Build genes dendrogram
dst <- dist(t(scale(t(Smooth.curve.matrix))), method = "manhattan")
hc <- hclust(dst, method = "ward.D") #"ward.D"

# Seriate the dendrogram
dend <- as.dendrogram(hc)
dend <- dendextend::seriate_dendrogram(dend, dst, method="GW")

# Find 15 cluster based on hierarchical clustering
clusters <- cutree(hc, k=10)

Gene.dynamique <- data.frame(Gene= names(clusters),
                             pval= DV.Axis.genes.FDR.filtered$pval,
                             qval=DV.Axis.genes.FDR.filtered$qval,
                             num_cells_expressed=DV.Axis.genes.FDR.filtered$num_cells_expressed,
                             Waves = as.numeric(clusters)) %>% arrange(Waves)

row.names(Gene.dynamique) <- Gene.dynamique$Gene
Gene.dynamique$Gene.Clusters <- paste0("Clust.",Gene.dynamique$Waves)

write.table(Gene.dynamique, "./Progenitors/Gene.dynamique.csv", sep = ";", quote = F)

Plot gene/cell heatmap

anno.colors <- list(Domaines = c(Clust.1="#83c3b8", Clust.2="#009fda", Clust.3="#3e69ac", Clust.4="#e46b6b", Clust.5="#e3c148", Clust.6="#b7d174", Clust.7="#68b041"),
                    Gene.Clusters = c(Clust.1 ="#ec756d" , Clust.2="#ebcb2e", Clust.3="#b79f0b", Clust.4="#7293c8", Clust.5="#3ca73f",
                                      Clust.6="#9ec22f", Clust.7="#cc3a1b", Clust.8="#4cabdc", Clust.9="#cc8778", Clust.10="#4784a2"))

pheatmap::pheatmap(Smooth.curve.matrix,
                   cluster_rows= as.hclust(dend),
                   scale = "row",
                   cluster_cols = F,
                   gaps_col = cumsum(as.numeric(table(Domaines.Clust$Domaines))),
                   annotation_row = Gene.dynamique %>% dplyr::select(Gene.Clusters),
                   annotation_col = Domaines.Clust,
                   annotation_colors = anno.colors,
                   show_colnames = F,
                   show_rownames = T,
                   fontsize_row = 2,
                   color = rev(brewer.pal(11,"RdBu")),
                   breaks = seq(-2.5,2.5, length.out = 11),
                   main = "Genes expression along Dorso-Ventral axis")

Manuscript Fig. 5C

##Plot gene clusters trends

source("./functions/GenesTrendPlots.R")
Clusters.trend(AP.data,
               Which.cluster = 1:10,
               clust.list = clusters,
               group.by = "global",
               span = 1,
               Smooth.method = "auto",
               Use.scale.data = T)

Manuscript Fig. 5C

Assign domain identity

We assign domain identity based on clusters’ transcriptional profile by setting boundaries over pseudo-dv score

# Set the boundary over speudotime score
new_data$cluster <- Domaines.Clust$Domaines
Infered.Domain.boundary <- aggregate(DorsoVentral.Score ~ cluster, new_data, max) %>% pull(DorsoVentral.Score)

# Assign identity based on the position of the cell on the pseudo-dv axis
Domaine.Ident <- sapply(AP.data@meta.data$DorsoVentral.Score,
                        function(x){ if(x<Infered.Domain.boundary[1]){ x = "Sub.Pallium.1"
                        } else if(x> Infered.Domain.boundary[1] & x< Infered.Domain.boundary[2]){ x ="Sub.Pallium.2"
                        } else if(x> Infered.Domain.boundary[2] & x< Infered.Domain.boundary[3]){ x = "Sub.Pallium.3" 
                        } else if(x> Infered.Domain.boundary[3] & x< Infered.Domain.boundary[4]){ x = "Ventral.Pallium"
                        } else if(x> Infered.Domain.boundary[4] & x< Infered.Domain.boundary[5]){ x = "lateral.Pallium.1"
                        } else if(x> Infered.Domain.boundary[6]){ x = "Dorsal.Pallium"
                        } else x="lateral.Pallium.2"})

# Transfert the identity to the Seurat object
AP.data@meta.data$Domaine <- Domaine.Ident
AP.data <- SetAllIdent(AP.data, id = "Domaine")

Plot.Genes.trend(AP.data,
                 genes = c("Gsx2", "Dbx1", "Gm29260", "Tfap2c", "Emx1", "Lrrn1"),
                 Use.scale.data = F)

Manuscript Fig. 5D

DimPlot(AP.data,
         group.by = "Domaine",
         reduction.use = "spring",
         dim.1 = 1,
         dim.2 = 2,
         do.label=F,
         label.size = 4,
         no.legend = F,
         cols.use = tolower(c("#68B041", "#E3C148", "#B7D174", "#83C3B8", "#009FDA", "#3E69AC", "#E46B6B"))
         )

Plot Figure 5B

#load full dataset
Allcells.data <- readRDS("./QC.filtered.cells.RDS")

#Transfer the identities
Rename.Clust <-  function(Clustdata, RawQCdata) {
  unClustered.cells <- RawQCdata@meta.data$Barcodes
  RawQCdata <- SetIdent(RawQCdata, cells.use = unClustered.cells, ident.use = "All.Unclustered.Cells")
  
  for(i in unique(Clustdata@meta.data$Domaine)){
  New.ident <- i
  Barcodes <- rownames(subset(Clustdata@meta.data, Clustdata@meta.data$Domaine == i))
  print(paste0("Cluster_",i,": ",length(Barcodes), " Cells"))
  Barcodes <- Barcodes[Barcodes %in% rownames(RawQCdata@meta.data)]
  RawQCdata <- SetIdent(RawQCdata, cells.use = Barcodes ,ident.use = paste0("AP.",i))
  }
  return(RawQCdata)
}

Allcells.data <- Rename.Clust(Clustdata = AP.data, RawQCdata = Allcells.data)

## [1] "Cluster_Sub.Pallium.2: 363 Cells"
## [1] "Cluster_Ventral.Pallium: 312 Cells"
## [1] "Cluster_Sub.Pallium.3: 197 Cells"
## [1] "Cluster_lateral.Pallium.1: 219 Cells"
## [1] "Cluster_Dorsal.Pallium: 201 Cells"
## [1] "Cluster_lateral.Pallium.2: 113 Cells"
## [1] "Cluster_Sub.Pallium.1: 241 Cells"

colors <-  c("#969696",tolower(c("#68B041", "#E3C148", "#B7D174", "#83C3B8", "#009FDA", "#3E69AC", "#E46B6B")))

DimPlot(Allcells.data,
        reduction.use = "spring", 
        dim.1 = 1,
        dim.2 = 2,
        do.label=T,
        label.size = 2,
        no.legend = T,
        cols.use = colors)

Manuscript Fig. 5B

rm(list = ls()[!ls() %in% c("AP.data")])

Plot the representative gene of the Figure S6

# Load custom plotting functions
source("./functions/GenesTrendPlots.R")

Cluster 6

Plot.Genes.trend(AP.data,
                 genes = c("Dlx1", "Ano1", "Dlx2", "Olig2"),
                 Use.scale.data = F)

Manuscript Fig.S6A

Cluster 2

Plot.Genes.trend(AP.data,
                 genes = c("Otx2", "Zbtb20", "Six3", "Ascl1"),
                 Use.scale.data = F)

Manuscript Fig.S6B

Cluster 3

Plot.Genes.trend(AP.data,
                 genes = c("Meis1", "Rorb", "Epha3"),
                 Use.scale.data = F)

Manuscript Fig.S6C

Cluster 1

Plot.Genes.trend(AP.data,
                 genes = c("Sfrp2", "Eva1c", "Sema5a", "Fat4"),
                 Use.scale.data = F)

Manuscript Fig.S6D

Cluster 9

Plot.Genes.trend(AP.data,
                 genes = c("Dmrta2","Lypd6", "Pax6", "Dmrta1"),
                 Use.scale.data = F)

Manuscript Fig.S6E

Cluster 5 ?

Plot.Genes.trend(AP.data,
                 genes = c("Emx1", "Arx", "Lmo3"),
                 Use.scale.data = F)

Manuscript Fig.S6F

Cluster 9

Plot.Genes.trend(AP.data,
                 genes = c("Fezf2", "Emx2", "Sp8"),
                 Use.scale.data = F)

Manuscript Fig.S6G

Transfert these ident on the full dataset (Manuscript Fig. 2A)

#load full dataset
Allcells.data <- readRDS("./Clustered.cells.RDS")

#Transfer the identities
Rename.Clust <-  function(Clustdata, RawQCdata) {
  for(i in unique(Clustdata@meta.data$Domaine)){
    New.ident <- i
    Barcodes <- rownames(subset(Clustdata@meta.data, Clustdata@meta.data$Domaine == i))
    print(paste0("Cluster_",i,": ",length(Barcodes), " Cells"))
    Barcodes <- Barcodes[Barcodes %in% rownames(RawQCdata@meta.data)]
    RawQCdata <- SetIdent(RawQCdata, cells.use = Barcodes ,ident.use = paste0("AP.",i))
  }
  return(RawQCdata)
}

Allcells.data <- Rename.Clust(Clustdata = AP.data, RawQCdata = Allcells.data)

## [1] "Cluster_Sub.Pallium.2: 363 Cells"
## [1] "Cluster_Ventral.Pallium: 312 Cells"
## [1] "Cluster_Sub.Pallium.3: 197 Cells"
## [1] "Cluster_lateral.Pallium.1: 219 Cells"
## [1] "Cluster_Dorsal.Pallium: 201 Cells"
## [1] "Cluster_lateral.Pallium.2: 113 Cells"
## [1] "Cluster_Sub.Pallium.1: 241 Cells"

colors2 <-  c("#969696",
              tolower(c("#68B041", "#E3C148", "#B7D174", "#83C3B8", "#009FDA", "#3E69AC", "#E46B6B")),
              "#ec756d", "#c773a7", "#7293c8", "#b79f0b", "#3ca73f","#31b6bd",
              "#ebcb2e", "#9ec22f", "#a9961b", "#cc3a1b", "#cc8778" , "#d14c8d", "#4cabdc", "#5ab793", "#e7823a","#e6bb9b", "#046c9a", "#4784a2" , "#4990c9")

DimPlot(Allcells.data,
        reduction.use = "spring", 
        dim.1 = 1,
        dim.2 = 2,
        do.label=T,
        label.size = 2,
        no.legend = T,
        cols.use = colors2)

Manuscript Fig. 2A

Allcells.data <- StashIdent(Allcells.data, save.name = "Cluster.ident")
saveRDS(Allcells.data, "./Clustered.cells.RDS")
saveRDS(AP.data, "./AP.data.RDS")

rm(list = ls()[!ls() %in% c("AP.data")])

Perform iterative graph-based clustering to find if discrete pallial sub-domain exist

Extract pallial progenitors

# Filter subpallial progenitors
AP.ident <- unique(AP.data@ident)
AP.data <-  SubsetData(AP.data,
                       ident.use = grep("Sub", AP.ident, value = T, invert = T),
                       subset.raw = T,
                       do.clean = F)

# Reset the pseudo-DV score to [0,1]
score <- AP.data@meta.data$DorsoVentral.Score
AP.data@meta.data$DorsoVentral.Score <- (score - min(score)) / (max(score) - min(score))

Perform iterative clustering implemented in scrattch.hicat package

For more detail on the scrattch.hicat please refer to the package page.

Prepare the dataset for clustering with scrattch.hicat

Independent gene filtering

We decide to exclude Cell cycle, ribosomal and mitochondrial associated genes, as well as Xist for the clustering step. We also retained gene detected for which transcripts have been detected in at least 3 cells

# Exclude cell cycle associated genes
CCgenes <- as.character(read.table("./Progenitors/CellCycleGenes.csv", sep = "\t", header = F)[,1])
#Exclude genes detected in less than 3 cells
num.cells <- Matrix::rowSums(AP.data@data > 0)
genes.use <- names(x = num.cells[which(x = num.cells >= 3)])

GenesToRemove <- c(grep(pattern = "(^Rpl|^Rps|^Mrp)", x = genes.use, value = TRUE),
                   grep(pattern = "^mt-", x = genes.use, value = TRUE),
                   "Xist", CCgenes)

genes.use <- genes.use[!genes.use %in% GenesToRemove]

Normalization

dgeMatrix_count <- as.matrix(AP.data@raw.data)[rownames(AP.data@raw.data) %in% genes.use,]
dgeMatrix_cpm <- cpm(dgeMatrix_count)
norm.dat <- log2(dgeMatrix_cpm + 1)

norm.dat <- Matrix(norm.dat, sparse = TRUE)
Data.matrix <- list(raw.dat=dgeMatrix_count, norm.dat=norm.dat)
attach(Data.matrix)

rm(list = ls()[!ls() %in% c("AP.data", "Data.matrix")])

Select feature used to remove PCs which would correlate with unwanted source of variations

#Dimension Filtering
gene.counts <- log2(colSums(as.matrix(Data.matrix$norm.dat) > 0))
nUMI <- log2(colSums(Data.matrix$raw.dat))
perctMito <- AP.data@meta.data$percent.mito
perctRibo <- AP.data@meta.data$percent.ribo
PseudoDVscore <- AP.data@meta.data$DorsoVentral.Score

rm.eigen <- as.matrix(cbind(gene.counts,
                            nUMI,
                            perctMito,
                            perctRibo,
                            PseudoDVscore))

row.names(rm.eigen) <- names(gene.counts)

colnames(rm.eigen) <- c("log2nGenes",
                        "log2nUMI",
                        "perctMito",
                        "perctRibo",
                        "PseudoDVscore")

rm(gene.counts, nUMI, perctMito, perctRibo,PseudoDVscore)

Iterative clustering

# Setting Parameters for iterative clustering
de.param <- de_param(padj.th     = 0.01, 
                     lfc.th      = 0.9,
                     low.th      = 1, 
                     q1.th       = 0.25, 
                     q2.th       = NULL,
                     q.diff.th   = 0.7,
                     de.score.th = 100,
                     min.cells = 10)

Run the iterative clustering

The default iter_clust function use in this version of the scrattch.hicat package does not allow to set the k.param argument. We modify this function to allow this argument to be set to other values.

# Perform the iterative clustering 
iter.result <- iter_clust(norm.dat, 
                          counts = raw.dat,
                          dim.method = "pca",
                          max.dim = 15,
                          de.param = de.param,
                          type = "undirectional",
                          rm.eigen = rm.eigen,
                          k.param = 15,
                          rm.th = 0.7,
                          vg.padj.th = 0.5,
                          method = "louvain",
                          prefix = "test-iter_clust",
                          verbose = F)

## [1] "test-iter_clust"
##   Finding nearest neighbors...DONE ~ 0.026 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.048 s
##   Build undirected graph from the weighted links...DONE ~ 0.057 s
##   Run louvain clustering on the graph ...DONE ~ 0.034 s
##   Return a community class
##   -Modularity value: 0.6822028 
##   -Number of clusters: 12[1] "test-iter_clust.1"
##   Finding nearest neighbors...DONE ~ 0.015 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.044 s
##   Build undirected graph from the weighted links...DONE ~ 0.046 s
##   Run louvain clustering on the graph ...DONE ~ 0.021 s
##   Return a community class
##   -Modularity value: 0.6939065 
##   -Number of clusters: 12[1] "test-iter_clust.2"
##   Finding nearest neighbors...DONE ~ 0 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.008 s
##   Build undirected graph from the weighted links...DONE ~ 0.009 s
##   Run louvain clustering on the graph ...DONE ~ 0.002 s
##   Return a community class
##   -Modularity value: 0.7585217 
##   -Number of clusters: 7

Merge clusters not seperable by DEGs

# Merge clusters which are not seperable by DEGs because they would have been separated at a early step of the process
rd.dat <- t(norm.dat[iter.result$markers,])
merge.result <- merge_cl(norm.dat, 
                         cl = iter.result$cl, 
                         rd.dat = rd.dat,
                         de.param = de.param)

cat(length(unique(merge.result$cl))," Clusters\n")

## 2  Clusters

cat(length(unique(merge.result$markers)), " DE genes")

## 10  DE genes

Import cell labels into the Seurat obj

AP.data@ident <- as.factor(merge.result$cl)
AP.data@meta.data$Pallial.subDomaine <-as.character(AP.data@ident)

colors <-  c("#969696", "#e7823a")

DimPlot(AP.data,
        reduction.use = "spring",
        dim.1 = 1,
        dim.2 = 2,
        do.label=T,
        pt.size = 2,
        label.size = 4,
        no.legend = F,
        cols.use = colors)

Manuscript element Fig. 2C

Find sub-domain marker genes

# Find genes which best distinguish the subdomain
Allmarkers <- FindAllMarkers(object = AP.data,
                             test.use = "roc",
                             min.pct = 0.3,
                             logfc.threshold = 0.6,
                             print.bar = T,
                             only.pos = T)

Allmarkers$q.diff <- abs(Allmarkers$pct.1 - Allmarkers$pct.2)/max(Allmarkers$pct.1, Allmarkers$pct.2)
Allmarkers$Specificity.index <- Allmarkers$avg_logFC/ (1-Allmarkers$q.diff)

Selected.markers <- Allmarkers %>%
  filter(Specificity.index > 1) %>% pull(gene)

# Prepare annotation for hicat pipeline
colors <-  c("#969696", "#e7823a")
colorsident <- cbind(ident = unique(as.character(AP.data@ident)),
                     colors = colors,
                     id = unique(as.character(AP.data@ident)))

#Create annotation data.frame
anno.df <- as.data.frame(cbind(
  sample_name = row.names(AP.data@meta.data),
  primary_type_id = colorsident[match(as.character(AP.data@ident), colorsident[,1]),3],
  primary_type_label = as.character(AP.data@ident),
  primary_type_color = colorsident[match(as.character(AP.data@ident), colorsident[,1]),2]
))


#Plot cells barplot
data <- cbind(sample_name = colnames(AP.data@data),
              as.data.frame(t(as.matrix(AP.data@data))))

sample_bar_plot(data, 
                anno.df, 
                genes = Selected.markers,
                group_order = c("1","3"),
                grouping = "primary_type",
                log_scale = FALSE,
                font_size = 7,
                label_height = 10,
                label_type = "angle",
                bg_color ="#f7f7f7")

Session Info

#date
format(Sys.time(), "%d %B, %Y, %H,%M")

## [1] "09 novembre, 2020, 10,20"

#Packages used
sessionInfo()

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
##  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
##  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] splines   stats4    parallel  stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] limma_3.42.0          Rphenograph_0.99.1    igraph_1.2.5         
##  [4] matrixStats_0.55.0    viridis_0.5.1         viridisLite_0.3.0    
##  [7] wesanderson_0.3.6     RColorBrewer_1.1-2    gridExtra_2.3        
## [10] reshape_0.8.8         ggExtra_0.9           scrattch.vis_0.0.210 
## [13] purrr_0.3.3           ggbeeswarm_0.6.0      dplyr_0.8.3          
## [16] scrattch.hicat_0.0.16 cluster_2.1.0         princurve_2.1.4      
## [19] monocle_2.14.0        DDRTree_0.1.5         irlba_2.3.3          
## [22] VGAM_1.1-2            Biobase_2.46.0        BiocGenerics_0.32.0  
## [25] Seurat_2.3.4          Matrix_1.2-17         cowplot_1.0.0        
## [28] ggplot2_3.2.1        
## 
## loaded via a namespace (and not attached):
##   [1] snow_0.4-3           backports_1.1.5      Hmisc_4.3-0         
##   [4] plyr_1.8.4           lazyeval_0.2.2       densityClust_0.3    
##   [7] fastICA_1.2-2        digest_0.6.25        foreach_1.4.7       
##  [10] htmltools_0.5.0      lars_1.2             gdata_2.18.0        
##  [13] magrittr_1.5         checkmate_1.9.4      gclus_1.3.2         
##  [16] mixtools_1.1.0       ROCR_1.0-7           R.utils_2.9.0       
##  [19] docopt_0.6.1         colorspace_1.4-1     ggrepel_0.8.1       
##  [22] xfun_0.18            sparsesvd_0.2        crayon_1.3.4        
##  [25] jsonlite_1.7.0       zeallot_0.1.0        survival_2.44-1.1   
##  [28] zoo_1.8-6            iterators_1.0.12     ape_5.3             
##  [31] glue_1.4.1           registry_0.5-1       gtable_0.3.0        
##  [34] kernlab_0.9-29       prabclus_2.3-1       DEoptimR_1.0-8      
##  [37] scales_1.1.0         pheatmap_1.0.12      bibtex_0.4.2        
##  [40] miniUI_0.1.1.1       Rcpp_1.0.5           metap_1.1           
##  [43] dtw_1.21-3           xtable_1.8-4         htmlTable_1.13.2    
##  [46] reticulate_1.13      foreign_0.8-72       bit_4.0.4           
##  [49] proxy_0.4-23         mclust_5.4.5         SDMTools_1.1-221.1  
##  [52] Formula_1.2-3        tsne_0.1-3           htmlwidgets_1.5.1   
##  [55] httr_1.4.1           FNN_1.1.3            gplots_3.0.1.1      
##  [58] fpc_2.2-3            acepack_1.4.1        modeltools_0.2-22   
##  [61] ica_1.0-2            farver_2.0.1         pkgconfig_2.0.3     
##  [64] R.methodsS3_1.7.1    flexmix_2.3-15       nnet_7.3-14         
##  [67] labeling_0.3         later_1.0.0          tidyselect_0.2.5    
##  [70] rlang_0.4.7          reshape2_1.4.3       munsell_0.5.0       
##  [73] tools_3.6.3          ggridges_0.5.1       fastmap_1.0.1       
##  [76] evaluate_0.14        stringr_1.4.0        yaml_2.2.1          
##  [79] npsurv_0.4-0         knitr_1.26           bit64_4.0.2         
##  [82] fitdistrplus_1.0-14  robustbase_0.93-5    caTools_1.17.1.2    
##  [85] RANN_2.6.1           dendextend_1.12.0    pbapply_1.4-2       
##  [88] nlme_3.1-141         mime_0.7             slam_0.1-46         
##  [91] R.oo_1.23.0          hdf5r_1.3.2.9000     compiler_3.6.3      
##  [94] rstudioapi_0.11      beeswarm_0.2.3       png_0.1-7           
##  [97] lsei_1.2-0           tibble_2.1.3         stringi_1.4.6       
## [100] highr_0.8            lattice_0.20-41      HSMMSingleCell_1.6.0
## [103] vctrs_0.2.0          pillar_1.4.2         lifecycle_0.1.0     
## [106] combinat_0.0-8       Rdpack_0.11-0        lmtest_0.9-37       
## [109] data.table_1.12.6    bitops_1.0-6         seriation_1.2-9     
## [112] gbRd_0.4-11          httpuv_1.5.2         R6_2.4.1            
## [115] latticeExtra_0.6-28  TSP_1.1-10           promises_1.1.0      
## [118] KernSmooth_2.23-15   vipor_0.4.5          codetools_0.2-16    
## [121] MASS_7.3-53          gtools_3.8.1         assertthat_0.2.1    
## [124] withr_2.1.2          qlcMatrix_0.9.7      mgcv_1.8-33         
## [127] diptest_0.75-7       doSNOW_1.0.18        grid_3.6.3          
## [130] rpart_4.1-15         tidyr_1.0.0          class_7.3-17        
## [133] rmarkdown_2.5        segmented_1.0-0      Rtsne_0.15          
## [136] shiny_1.4.0          base64enc_0.1-3

---

1. Institute of Psychiatry and Neuroscience of Paris, INSERM U1266, 75014, Paris, France↩